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Measures of multiple spike train synchrony are essential in order to study issues such as spike 
timing reliability, network synchronization, and neuronal coding. These measures can broadly be 
divided in multivariate measures and averages over bivariate measures. One of the most recent 
bivariate approaches, the ISI-distance, employs the ratio of instantaneous interspike intervals (ISIs). 
In this study we propose two extensions of the ISI-distance, the straightforward averaged bivariate 
ISI-distance and the multivariate ISI-diversity based on the coefficient of variation. Like the orig- 
inal measure these extensions combine many properties desirable in applications to real data. In 
particular, they are parameter free, time scale independent, and easy to visualize in a time-resolved 
manner, as we illustrate with m vitro recordings from a cortical neuron. Using a simulated network 
of Hindemarsh-Rose neurons as a controlled configuration we compare the performance of our meth- 
ods in distinguishing different levels of multi-neuron spike train synchrony to the performance of six 
other previously published measures. We show and explain why the averaged bivariate measures 
perform better than the multivariate ones and why the multivariate ISI-diversity is the best per- 
former among the multivariate methods. Finally, in a comparison against standard methods that 
rely on moving window estimates, we use single-unit monkey data to demonstrate the advantages 
of the instantaneous nature of our methods. 

Keywords: time series analysis; spike trains; reliability; clustering; neuronal coding; neuronal networks; 
synchronization; ISI-distance 



INTRODUCTION 



Estimating the degree of synchrony within a set of 
spike trains is a common task within two different sce- 
narios. In the first scenario spike trains are recorded 
successively from only one neuron. The most common 
application is to quantify the reliability of the neuronal 
response upon repeated presentations of the same stim- 
ulus fvA [26| . In the second scenario spike trains are 
derived from simultaneous recordings of a population of 
neurons 0, Q. A typical application for such data is 
the analysis of synchronization and desynchronization in 
large neural networks 

Measures designed to estimate single-neuron reliability 
and multi-neuron synchrony are either multivariate or bi- 
variate. Multivariate approaches evaluate all spike trains 
at the same time, e.g., by measuring the normalized vari- 
ance of pooled, exponentially convolved spike trains p^ . 
or by exploiting the deviation of pooled spike train statis- 
tics from the one obtained for a Poisson process [l^ . Bi- 
variate measures aim at the quantification of synchrony 
between just two spike trains, but can also be used in 
a multivariate context by defining the total synch rony 
as the average over all pairwise synchrony- values [l3|. 
Prominent examples are the cost-based distance intro- 
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duced in Victor and Purpura [28[ , the Euclidean distance 
proposed in van Rossur n 12711 . cross correlation of spike 
trains after filtering 0, |22{. and event synchronization 

Most recently, the bivariate ISI-distance has been in- 
troduced as a simple approach that extracts informa- 
tion from the interspike intervals (ISIs) by evaluating 
the ratio of instantaneous firing rates [ij]. The ISI- 
distance combines a variety of properties that make it 
ideally suited for applications to real data. In particu- 
lar, it serves as an excellent means to visualize the oc- 
currence of spiking patterns. In contrast to most of the 
other measures that rely on a parameter determining the 
time scale, the ISI-distance does not depend on any time 
scales and is thus free of parameters. In a comparison 
with previously published approaches on spike trains ex- 
tracted from a simulated Hindemarsh-Rose network, the 
ISI-distance matched the performance of the best time 
scale optimized measure [14 1. 

In this study we propose two extensions of the ISI- 
distance, the straightforward averaged bivariate ISI- 
distance and the multivariate ISI-diversity based on the 
coefficient of variation. Both of these extensions inherit 
the fundamental properties of the measure they are based 
on. In particular, they are parameter free, time scale in- 
dependent and easy to visualize in a time-resolved man- 
ner. We apply the two measures to in vitro recordings 
of a cortical cell from a rat and show that they serve as 
an excellent means to track spiking patterns in multiple 
spike trains. 
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In order to better characterize the two extensions of 
the ISI-distance, we consider them in the context of 
the groups of averaged bivariate and multivariate mea- 
sures to which they each belong. In particular, we com- 
pare their performance in distinguishing different lev- 
els of multi-neuron spike train synchrony to the perfor- 
mance of six other previously published averaged bivari- 
ate and multivariate measures. As a controlled configu- 
ration we use a simulated network of Hindemarsh-Rose 
neurons with predefined clustering [l^, [l9| . In this com- 
parison, the averaged bivariate measures perform better 
than the multivariate measures. Finally, we present a 
typical application using single-unit data from a visual 
attention task in a macaque monkey p^ . Comparing 
the ISI-measures to standard methods that rely on mov- 
ing window estimates, we stress the advantages of their 
instantaneous nature of our proposed methods. 

The remainder of the paper is organized as follows: In 
Section we present a more detailed description of the 
original bivariate ISI-distance (Section III Ap . the aver- 
aged bivariate ISI-distance fSection III Bp . and the mul- 
tivariate ISTdiversity (Section III Cp . Furthermore, we 
explain how to evaluate the significance of the values ob- 
tained in order to study properties of neuronal coding 
fSection III Pp . All of these measures are applied to in 
vitro recordings from a cortical neuron in Section [ml 
Subsequently, in Section HVl a performance comparison is 
carried out on a network of Hindemarsh-Rose neurons. 
Finally, in Section |V] a typical field data application is 
shown before the conclusions are drawn in Section IVII 
In the Appendix we derive analytical results for Poisson 
processes in Section |3 introduce the in vitro recordings 
of cortical cells, the simulated Hindemarsh-Rose time se- 
ries, and the in vivo monkey data in Sections IB II IB 21 
and IB 31 respectively, and describe the previously pub- 
lished approaches to measure spike train synchrony in 
Section O 



A. Bivariate ISI-distance 

To define a time-resolved, symmetric, and scale- 
invariant measure of the relative firing rate pattern 
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This quantity becomes for identical ISI in the two spike 
trains, and approaches —1 and 1, respectively, if the first 
(second) spike train is much faster than the other. 

Finally, the absolute ISI-distance is integrated over 
time: 



Dj 



dt\h2{t)V 
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Identical spike trains yield a value of zero if there is no 
phase lag between the spike trains (for ways to deal with 
non-zero phase lags please confer the discussion in Sec- 
tion |Vl|. Irrespective of any phase lag the value zero is 
also obtained for periodic spike trains with the same pe- 
riod (i.e., all ISIs of the two spike trains have the same 
length) . For the more general case of periodic spike trains 
with periods p and q (with p < q) a value of 



Di = l- p/q 



(5) 



is obtained. For the special case that p and q are integer 
values, this is the classical p : q synchronization. 



B. Averaged bivariate ISI-distance 



II. THE BIVARIATE ISI-DISTANCE AND ITS 
EXTENSIONS 

We denote with {tf} — ti,...,t^j^ the spike times 
and with Af„ the number of spikes for neuron n with 
n = 1,...,N. As a first step common to all the follow- 
ing instantaneous measures for each neuron the value of 
the current interspike interval is assigned to each time 
instant [sH 

x^Siit) = min(<r |tr >t)- ma^m7 < t) t^ < t < . 

(1) 

From the instantaneous ISIs the average instantaneous 
rate can be defined as 

1 ^ 1 



We again start by assigning the ISI for each spike train 
{<:"} with n — 1, ...,N as in Eq. [T] and proceed by cal- 
culating the instantaneous average A(t) over all pairwise 
absolute ISI-ratios |/m„(t)| (cf. Eq. [3]) 

^ N N 

^ ' n—1 m— n+1 

Averaging over time yields 

I?? r dtAit). (7) 

-'- Jt=o 

The same kind of time-resolved visualization as in the 
bivariate case is possible, because the two averages com- 
mute. Instead of averaging over pairwise ISI-distances, 
which themselves are obtained by averaging over time, we 
can first average over pairwise instantaneous ISI-ratios 
and then average over time. 
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C. Multivariate ISI-diversity based on coefficient of 
variation 

We derive a multivariate measure by taking the instan- 
taneous coefficient of variation taken across all neurons 
n = 1, iV at any given instant in time 



Cv{t) 



{t))n 



"ISI 



and again integrate over time: 



dtCv{t). 



(8) 



(9) 



The coefficient of variation is probably the most widely 
used measure of variability. Since it relates the standard 
deviation to the mean (and therefore is dimensionless), 
it allows for broader comparisons also between data sets 
with different units or different means. For identical spike 
trains and for periodic spike trains with the same period 
it obtains the same lower bound of zero as D/ and D° 
but unlike them it lacks an upper bound. 

The coefficient of variation can be very sensitive to out- 
liers in the instantaneous ISI distribution. Since this dis- 
tribution is bounded below by zero and typically skewed 
towards larger values, we first log-transform the data and 
then calculate the back-transformed coefficient of vari- 
ance according to Hopkins et al. 



C'°«(i) = (ln(xrsi(<))>„exp( 



<7{\n{x^^,{t))) 
(ln(.T«i(t)))„- 



(10) 



D. ISI-distances and neuronal coding 

The ISI-distance and its extensions provided here, as 
well as the other spike train distances, are not designed 
to directly address the question of what is coded by the 
neurons and in which way this is done, as for example 
is the case for many approaches that rely on spike trig- 
gering (23j . However, they allow examining properties 
such as the reliability of a neuron to the same stimulus 
or similarity across a population in a single presentation 
of a stimulus which can help to infer and characterize the 
neuronal response. For that purpose it is important to 
compare the values obtained to a baseline level indicating 
the results that would be obtained if different spike trains 
are assumed to be independent. One of the most com- 
mon reference models in spike train analysis is the Pois- 
son process, an example for a renewal process for which 
at each time instant the probability of a spike is indepen- 
dent of the occurrence of previous spikes. In Appendix 
El we investigate the case of multiple Poisson processes 
and derive analytical results for the ISI-distance D/, the 
averaged bivariate ISI-distance 1?° and the ISI-diversity 
DY all of which can be used as benchmarks for this type 
of random spike trains. 
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FIG. 1: (color online) Input and output spike trains recorded 
from a cortical cell. The same input patterns were delivered 
with three different firing rates that were each repeated five times 
(Trials 1-5: faster input, Trials 6-10: matched input, Trials 11- 
15: slower input). The input is shown in blue and the output in 
red. 



III. ILLUSTRATION OF THE ISI-MEASURES 
USING IN VITRO RECORDINGS FROM A 
CORTICAL CELL 



In an example application, the averaged bivariate ISI- 
distance -Df and the multivariate ISI-diversity Df are il- 
lustrated using in vitro whole-cell recordings taken from 
a cortical cell of a young Long- Evans rat. For a descrip- 
tion of the data see Appendix IB II Sets of synaptic in- 
puts were delivered to a spiking neuron using three dif- 
ferent input firing rates with five repetitions each. All 
recordings including the synaptic input trains are shown 
in Fig.[l]and their pairwise ISI-ratios (Eq. [3]) are depicted 
in Fig. H 

In the first five trials the input is too fast for the neuron 
to follow, while in the last five trials the neuron does not 
wait for the slower input and maintains its own, faster 
firing rate. In trials 6 — 10 the input firing rate matches 
the neuron's own rate, such that except for small de- 
viations 1 : 1 synchronization can be observed. The 
mean ISI-distance for the first five trials (first half only) 
is Di = 0.43 while it is D/ = 0.24 for the last five tri- 
als. However, the lowest distances are obtained for the 
matched inputs, for which all distances are smaller than 
0.1 with a mean distance of Di = 0.066. This resonance- 
like behavior is consistent with results recently reported 
in Haas et al. Q . 

Here we use the bivariate ISI-distance Dj to address 
the input-output relation for each trial; next we employ 
the averaged bivariate ISI-distance Dj and the multi- 
variate ISI-diversity I?™ (after log-transform, Eq. fTU]) to 
estimate the reliability of the output across repetitions 
of the same input, and the dependence on the input's 
temporal scaling. In Figs. [3] to [5] we investigate spiking 
patterns in the output spike trains of faster, temporally 
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FIG. 2: (color online) Pairwise input-output ISI-ratios for the 
fifteen trials shown in Fig. [T] The ISI-distances Di, defined as 
the average over the absolute values of the ISI-ratios (Eq. [4|, 
are shown on the right side of each trial. For the fast inputs 
(Trials 1-5) this is the ISI-distance averaged over the first half of 
the recording only, the part after the last input spike is ignored. 
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FIG. 4: (color online) Same as Fig. [3] but this time for the 
five output spike trains recorded during stimulation with speed- 
matched input (red traces 6— 10 in Fig.[l]). For this example the 
ISI-measures attain the low values = 0.09 and DT = 0.11. 
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FIG. 5: (color online) As Fig. [3] but this time for slower inputs. 
Here the ISI-measures are Df = 0.18 and = 0.19. 



FIG. 3: (color online) Relative firing patterns for output spike 
trains recorded using faster input (red traces 1 — 5 in Fig. [TJ. 
On top we depict the spike trains, followed by the instantaneous 
ISI-values according to Eq. [l] and their instantaneous average 
(black line). At the bottom we show the pairwise average (Eq. 
|6|) and the coefficient of variation (Eq. [8|). Dotted vertical lines 
mark beginning and end of the recording; dashed lines in the 
lower part indicate the average values over time, respectively. 
These are the ISI-measures which for this set of spike trains are 
equal to Dj = 0.17 and = 0.19, respectively. For the first 
half only both measures equal 0.13. 



matched and slower inputs, respectively. 

For the faster inputs (trials 1-5), which lasted only 
half the recording, the change in firing behavior at the 
transition can easily be tracked (Fig. In the begin- 
ning all output spike trains are fast and regular (cf. in- 
stantaneous ISI in Fig. [3]) and exhibit similar spiking 
patterns (cf. instantaneous variabilities A and Cy^ in 
Fig. [3]). Immediately after the last input spike, at ap- 



proximately 5 seconds, some of the output spike trains 
become very irregular. This results in high values of the 
instantaneous pairwise averages and the instantaneous 
coefficients of variation. After this transition period the 
variability of the spike trains returns to the initial low 
values, although now at a somewhat lower rate (higher 
instantaneous ISIs). 

Much lower ISI-measures are obtained for the matched 
input firing rates (trials 6-10, cf. Fig. [4]). Firing is more 
uniform and reliability is higher as reflected by the low 
values of the instantaneous bivariate averages and the in- 
stantaneous coefficients of variation. Occasionally these 
time-resolved measures reach 0, marking periods of per- 
fect global ISI synchronization, whereas in the second half 
they exhibit an increasing trend, which correctly reflects 
missing spikes in some of the spike trains. 

For slower input (trials 11-15) the ISI-measures yield 
again higher values similar to the ones for slower inputs 
(cf. Fig. [5]). However, for slower inputs the high vari- 
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ability is spread uniformly over the whole time interval, 
while the high distances for faster inputs are mostly due 
to the large irregularities at the termination of the input. 

In summary, this qualitative analysis with the ISI- 
measures can confirm previous observations on the de- 
pendence of reliability on input firing rate Q • Further- 
more, the time-resolved visualization used here allows 
for a more detailed analysis of the spiking patterns than 
measures that compress all the information into one final 
value can achieve. 



IV. COMPARISON OF MEASURES USING 
SIMULATED HINDEMARSH-ROSE TIME 
SERIES 

In this section we evaluate how well the ISI-measures 
and six previously published methods perform in distin- 
guishing different levels of multi-neuron spike train syn- 
chrony. In order to change the level of synchrony in a 
controlled manner, we constructed sets of spike trains 
from two clusters and gradually varied the relative con- 
tributions of these clusters. This configuration is reminis- 
cent of multi-unit recordings in which a given set of spike 
trains might represent different types of neurons or neu- 
rons responding to different tasks. In the first part of this 
section we show that the synchrony within such a non- 
homogeneous neuronal population is dependent on both 
the total number and the relative mixture of neurons. 
In the second part we define a criterion for comparing 
the performance of the various measures in distinguish- 
ing different sets of spike trains. 

The dataset consists of 26 spike trains taken from a 
network of Hindmarsh-Rose neurons. It is a subsample 
of the data previously used in Kreuz et al. ^14*1 . In a sim- 
ulation the network architecture was designed such that 
the 26 spike trains arise from two principal clusters with 
13 neurons each. This organization is reflected by the 
pairwise distance matrix shown in Fig. [5] for the bivari- 
ate ISI-distance D/. Note that the first cluster is more 
strongly connected than the second; the corresponding 
mean ISI-distances are (-D/)i = 0.17 and (-07)2 = 0.21, 
respectively. 

As a benchmark test, we compared the averaged bi- 
variate ISI-distance Dj (Eq. [7]) and the multivariate 
ISI-diversity (after log-transform, Eq. [TU)l with four 
averaged bivariate methods and two multivariate mea- 
sures. The averaged bivariate measures were comprised 
of the spike train distances Dy and Z?^ introduced by 
Victor and Purpura [2^ and van Rossum [2X], the simi- 
larity measure Dg proposed by Schreiber et al. [1^ and 
event synchronization Dq [21| . Multivariate measures in- 
cluded the reliabilities D'^ and Dip by Hunter et al. [l^l 
and Tiesinga [2^ , respectively. Four of these measures 
{Dy, D^, Dg and D^) rely on the choice of a parameter 
that sets the time scale. Following Kreuz et al. [3] these 
time scales were optimized with respect to the present 
task (see below). For a description of the underlying 
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FIG. 6: (color online) Pairwise distance matrix for the bivariate 
ISI-distance Dj applied to 26 time series from the Hindemarsh- 
Rose model. Neurons and their spike trains are labelled by '1' and 
'2' depending on their affiliation to the two clusters. Black lines 
separate the two intra-cluster and the inter-cluster submatrices. 
Their respective mean ISI-distances are shown. Note the higher 
mean ISI-distance Df of cluster 2 as compared to the one of 
cluster 1. 

model confer to Appendix IB 21 Details about the pre- 
viously published measures and their parameters can be 
found in Appendix [Cl 

To evaluate how well the measures are able to distin- 
guish different sets of spike trains, we performed 12 differ- 
ent runs with the total number of spike trains increasing 
from N = 2 to N = 13. In each run we varied the num- 
bers Ni and N2 of spike trains from the two clusters but 
kept the total number iV = iVi -|-iV2 of spike trains within 
the set constant. For any given N we started with spike 
trains from cluster 2 only. Then we successively replaced 
one spike train from this cluster by a spike train from 
cluster 1, until in the end the set consisted only of spike 
trains from cluster 1 . This set transition is parameterized 
by the set imbalance 



which varies from — 1 to 1 along the way. 

As a first example. Fig. [7] shows the two ISI-measures 
calculated for = 12 spike trains that are all taken from 
cluster 2 (set imbalance b — —1). Similar to the examples 
shown in the first part and to the visualization used for 
the bivariate ISI-distance, parts of the spike trains that 
are very similar or very different from each other are eas- 
ily identified. In the beginning, small ISI-measures reflect 
a short transient period of almost complete synchroniza- 
tion. Later, the ISI-measures highlight short-time ex- 
cursions that appear when transitions between short and 
long ISIs are temporally offset between different spike 
trains. 

A similar behavior can be observed for A^ = 12 spike 
trains from cluster 1 (set imbalance b = 1, Fig. [5]). In 
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FIG. 7: (color online) Averaged bivariate ISI-distance and mul- 
tivariate ISI-diversity for A'^ = 12 spike trains all from cluster 2 
corresponding to a set imbalance of 6 = —1. In this case the 
ISI-measures are _D? = 0.21 and = 0.36. 



FIG. 9: (color online) ISI-measures for the balanced case (fe = 
0) with 6 spike trains from cluster 1 (blue) and 6 spike trains 
from cluster 2 (red). Correspondingly, the ISI-measures are much 
higher: Df = 0.32 and = 0.51. 



Spike 
trains 



iSi 









il 'i ' i ' 'i'i 1 1 i' 'ii" "i" "i 

I 1 1 1 1 II III II 1 III 1 1' 

II 1 1 1 [1 III Ii 1 1 II II II nil 1 nil 1 























0.5 1 1.5 

Data points 



FIG. 8: (color online) Example of the ISI-measures for a set 
imbalance 6=1 with all A'^ = 12 spike trains from cluster 1. For 
this stronger connected cluster smaller ISI-measures — 0.17 
and D'P = 0.28 are obtained. 



this example, the ISI-measures are lower than those ob- 
served for cluster 2, which is due to the stronger coupling 
within cluster 1 (cf. Fig. [S]) . Finally, a set in which both 
clusters contribute an equal number of 6 spike trains (set 
imbalance 6 = 0, Fig. [5]) yields ISI-measures that are 
much higher than either of the two previous examples. 

Fig. [inishows the results from an application of the av- 
eraged bivariate ISI-distance to one complete group of set 
transitions. For each curve, corresponding to a constant 
total number of spike trains N, a monotonic increase to 
the maximum value of D° is followed by a monotonic de- 
crease. There is a clear asymmetry - the values for the 
minimum negative set imbalance are higher than those 
for the maximum positive set imbalance - again reflect- 
ing the fact that the first cluster is more strongly coupled 
than the second one. Different traces, corresponding to 
different numbers of neurons in the set, overlap consider- 
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FIG. 10; (color online) Averaged bivariate ISI-distance versus 
the set imbalance (cf Eq. Ill|) for randomly chosen set tran- 
sitions and decreasing total numbers of spike trains A'^. The 
enlarged markers all lie on the line for = 12 spike trains and 
correspond to the examples shown in Figs. I7I9I 



ably; only the sets with small numbers of neurons attain 
mostly higher values. 

The latter effect can be observed for all measures, and 
can easily be understood for the averaged bivariate mea- 
sures. For these measures the average pair wise distance 
only depends on the distance matrix; in fact, for a sufR- 
ciently high total number of spike trains A'^, it only de- 
pends on the mean distances within each of the two clus- 
ters and on the mean distance between the two clusters 
(averaging over the respective submatrices yields {D)i, 
{0)2 and {D)i2 \ compare to Fig. [6|). 

For each total number of spike trains and each set im- 
balance, the number of pairs Npi and Np2 within each 
cluster, the number of pairs Npi2 between the two clus- 
ters and the total number of pairs Np can be combined 
with the respective mean distances to yield the expected 
distance 



E{Dx) = 



Npi{D 



XI 



Np2{Dx)2 



Npl2{Dx)l2 



(12) 
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FIG. 11: (color online) Expected values for the averaged bivari- 
ate ISI-distance versus the set imbalance (cf. Eq. Ilip for dif- 
ferent total numbers of spike trains TV. Values were calculated 
using Eq. [T2]and the mean values of the submatrices shown in 

Fig. m 



with X representing any of the bivariate measures. For 
the averaged bivariate ISI-distance the expected values 
E{Df) are shown in Fig. [11] Note that for normaUzed set 
imbalances ±1 always the mean intra-cluster distances 
{D)i = 0.17 and {D)2 = 0.21 (cf. Fig. © are obtained 
independently of the total spike train number. 

For a given pairwise distance matrix with given intra- 
and inter-cluster distances, the value E{Dx) only de- 
pends on the three coefficients. For a small total number 
of spike trains there are fewer intra-cluster than inter- 
cluster pairs, thus the distances tend to be higher. How- 
ever, for a constant set imbalance the ratio Npi2/{Npi + 
Np2) of inter-cluster and intra-cluster coefficients de- 
creases as the total number of spike trains N increases. 
With approaching infinity the ratio converges to 1 
from above, and thus the asymptotic curve for the ex- 
pected value connects the marked points in Fig. I 111 given 
by {Dx)2, (^jf>i/4+(i?x)2/4+(i^x>i2/2 (the relative 
amounts of the three submatrices marked in Fig. [6] for an 
infinitivcly large matrix) and {Dx)i (from left to right). 
However, since in reality the two intra-cluster submatri- 
ces and the inter-cluster submatrix are non-uniform, the 
averaged bivariate distances for a given random set tran- 
sition also depend on the standard deviation of the intra- 
and inter-cluster distances. Although this kind of reason- 
ing is only directly applicable to the bivariate measures, 
the observed behavior for the multivariate measures is 
very similar. Thus, in general the synchrony within a 
small non-homogeneous neuronal population is depen- 
dent not only on the relative mixture of clusters, but 
also on the total number of neurons. This fact is very 
important for multi-unit recordings, where typically only 
a relatively small number of neurons can be isolated, and 
where the fraction and the identity of the neurons coding 
for the presented stimulus at that moment in time is not 
known. 

In order to numerically evaluate the performance in 
distinguishing different sets, we randomly selected r = 
500 set transitions for each N = 2, 13, and applied all 
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FIG. 12: Averaged bivariate ISI-distance D": Distributions for a 
constant total number of A'^ = 6 spike trains. The set transition 
takes place from top to bottom with the set imbalance increasing 
from b — —1 to 6 = 1 (see right side). Histograms for adjacent 
normalized set imbalances are connected by dashed lines, and the 
accompanying numbers represent the value of the Kolmogorov- 
Smirnov statistic for the respective pair of distributions. 



measures to all sets. In this way for each total number 
of spike trains and every set imbalance we obtained a 
distribution of 500 values. Then, for every set transition 
we looked at the distributions of measure values for ad- 
jacent set imbalances only, and quantified their pairwise 
difference using Kolmogorov-Smirnov discrimination val- 
ues H. 

An example set transition for the averaged bivariate 
ISI-distance Dj with a constant total number of iV = 6 
spike trains is shown in Fig. [121 In this case the 
first and the last two pairs of distributions are disjunct, 
yielding Kolmogorov-Smirnov statistics equal to if = 1, 
whereas the intermediate distributions overlapped and 
correspondingly the Kolmogorov-Smirnov statistics at- 
tain values smaller than one. 

In Figs. [13] and [TJ] an overview of the Kolmogorov- 
Smirnov statistic for all constant total numbers of spike 
trains is shown for the averaged bivariate ISI-distance 
Dj and the Tiesinga dissimilarity D™ , respectively. Most 
curves exhibit a perfect distinction (maximum value K = 
1) for very large absolute imbalances, whereas the distri- 
butions for rather balanced sets of spike trains overlap 
considerably (smaller values of K). In all cases, mini- 
mum values are obtained for slightly negative set imbal- 
ances, again reflecting the different intra-cluster coupling 
strengths. However, the range of set imbalances for which 
the test discriminates perfectly is much broader for the 
averaged bivariate ISI-distance. 

From the ensemble of all of these traces, we defined 
our measure of performance, the set separation Ss, as the 
fraction of pairs of distributions that can be perfectly dis- 
tinguished, i.e., that yield a Kolmogorov-Smirnov statis- 
tic of A' = 1 (similar results are obtained for the 
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FIG. 13: (color online) Kolmogorov-Smirnov statistic K for the 
averaged bivariate ISI-distance D° and different total numbers 
of spike trains. The thick line corresponds to the distributions 
shown in Fig. 1121 Values are plotted against the mean set 
imbalance of the respective two distributions. For this measure 
the set separation is Ss ~ 0.54. 
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FIG. 15: Comparison of measures: Set separation Ss defined as 
the fraction of tests with a Kolmogorov-Smirnov statistic of K — 
1. The dashed line separates the averaged bivariate measures 
from the multivariate measures. 



V. APPLICATION TO MONKEY DATA FROM 
A VISUAL ATTENTION TASK 
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FIG. 14: (color online) Same as Fig. [13] but this time for the 
Tiesinga dissimilarity which yields a set separation of Ss = 
0.27. 



mean distinction, the average value over all Kolmogorov- 
Smirnov statistics). In order to assess the standard error 
we repeated the analysis five times using 100 realizations 
each time. In Fig. [15] we show the set separation Ss and 
its standard deviation for all measures. The most fun- 
damental result is that the averaged bivariate measures 
are better at distinguishing different set imbalances than 
the multivariate measures. The best distinctions are ob- 
tained for the averaged Victor-Purpura distance Dy and 
the van Rossum distance Z?^, which were each optimized 
with respect to the time scale parameter, followed by the 
averaged bivariate ISI-distance Dj, for which no opti- 
mization was needed. The multivariate ISI-diversity Z?™ 
after log-transform and without any optimization per- 
forms best among the multivariate measures. 



In the previous section we compared the time scale 
independent ISI-measures only against methods that ei- 
ther employ a parameter determining the time scale or 
are time scale independent themselves. However, some of 
the standard methods used in the literature rely on di- 
viding spike trains into disjunct or overlapping segments, 
which are then analyzed using a moving window analy- 
sis. In this last section we compare the ISI-measures with 
such an analysis using single-unit data from one neuron 
in area V4 of a rhesus macaque monkey during a visual 
attention task. The recordings and the task are described 
in Appendix IB 31 for more details see Mitchell et al. p^ . 
In this study Mitchell et al. calculated the firing rate 
and the Fano factor (a measure of response variability: 
the ratio of spike count variance to mean spike count) 
of one neuron for different attention conditions. The 
firing rate was smoothed using a Gaussian kernel with 
a = 12.5tos, whereas the Fano factor was calculated in 
non-overlapping 100 ms bins and then averaged across 
all visually responsive neurons for attended and unat- 
tended stimuli separately. As these two measures show 
(Fig. [16)) . attention leads to an increase in firing rate 
and a decrease in variability (increase in reliability) dur- 
ing the pause period in which the attended stimulus is 
positioned within the neuron's receptive field. 

In Fig. [16] we also depict the instantaneous firing rate 
R{t) (Eq. [5]) as well as the instantaneous estimates A{t) 
(Eq. ^ and Cy^t) (Eq. [TO]) for the averaged bivariate 
ISI-distance and the multivariate ISI-diversity, respec- 
tively. These estimates have the highest possible resolu- 
tion since they change with every occurrence of a single 
spike in one of the spike trains. However, in order to 
capture long term trends curves can be smoothed using 
a moving average filter with the desired order resulting in 
quantities R*{t), A*{t), and Cy'^*{t), respectively. The 
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(a parameter to be fixed beforehand), our proposed mea- 
sures are time-scale-independent and parameter-free. 

Furthermore, the method used above (fifth panel of 
Fig. [in]), the PSTH and other estimates of firing rate 
based on the summed spike count are invariant to shuf- 
fiing spikes among the spike trains. They yield the same 
high value for a window regardless of whether the spikes 
are spread across the different spike trains or whether 
all the spikes originate from the same spike train. In 
contrast, our instantaneous estimates take into account 
specific information from the individual spike trains by 
averaging over the intervals between their previous and 
their next spike. This method of averaging, in which the 
relevant unit of the average is the spike train and not the 
single spike, is closer to the actual notion of mean firing 
rate. 

Finally, for the case of Gaussian variability of the rate 
across trials it can be shown that the Fano factor is af- 
fected by changes in the mean rate ^] . In contrast, the in- 
stantaneous estimates of the ISI-distances are time-scale 
adaptive and thus do not depend on the mean rate. 



Time [s] 



VI. 



DISCUSSION 



FIG. 16: (color online) Single-unit monkey data from one neu- 
ron in area V4 during a visual attention task. In the two top 
panels we depict 96 spike trains and their instantaneous ISI, 
48 attended (red) and 48 unattended (blue) trials. The next 
three panels show firing rate estimates: the instantaneous fir- 
ing rate R defined as the instantaneous average over the inverse 
ISIs of the individual neurons, its moving average R* as well as 
the smoothed mean firing rate (denoted as "Rate"). The five 
bottom panels depict variability estimates based on the instan- 
taneous bivariate average A and the multivariate coefficient of 
variation Cy^ , their moving averages A* and C[?^*, as well as 
the Fano factor (denoted as "Fano"). In the latter panel a hori- 
zontal dashed line marks 1. Vertical black lines designate the be- 
ginning and the end of the interval in which the stimulus paused 
within the receptive field; the response starts approximately 400 
ms earlier when the stimulus enters the receptive field. 



moving average was applied in a manner that adapts to 
the local firing rate, i.e., the averaging was performed 
not over fixed time intervals but over fixed numbers of 
ISIs of the pooled spike train. The order was set to 50 
in order to be sensitive to the same time scales that were 
highlighted by Mitchell et al. 

A qualitative comparison between the standard meth- 
ods and our proposed measures (in particular, the mov- 
ing averages) reveals that the trends shown seem to be 
very similar. However, the instantaneous estimates are 
sensitive to additional features in the spike trains, since 
they exhibit maximum temporal resolution. Note also 
that the two kinds of approaches rely on different infor- 
mation. While the standard methods such as the Peri- 
Stimulus-Time-Histogram (PSTH) and the Fano factor 
capture features contained in windows of a certain size 



In this study we presented two extensions of the ISI- 
distance that measure the dissimilarity within a set of 
spike trains, which is typically derived either from suc- 
cessive recordings from only one neuron (e.g., upon re- 
peated presentations of the same stimulus), or from si- 
multaneous recordings of a population of neurons. The 
first extension employs the instantaneous average over 
the bivariate ISTratios whereas the second extension is 
multivariate and relies on the instantaneous coefficient of 
variation. Both of these extensions inherit the basic prop- 
erties of the ISI-distance: they are parameter free, time 
scale independent and easy to visualize in a time-resolved 
manner. Furthermore, they are conceptually simple and 
computationally fast. They also yield good results even 
on a rather limited number of short spike trains. In the 
first example we dealt with only 5 spike trains of about 
60 spikes each. 

As illustrated in the first part of this study in which 
we analyzed in vitro recordings of a cortical cell from 
a rat, the ISI-measures serve as an excellent means to 
track the occurrence of firing patterns in spike trains. 
Deviations from synchrony can be localized in time, thus 
rendering these methods a good choice for time-resolved 
applications on non-stationary neuronal dynamics. In 
the context of reliability estimations, this property allows 
the analyst to relate intervals of high or low synchrony 
to local stimulus features. 

In the second part of this study we showed that for 
small non-homogenous neuronal populations, the values 
of the averaged bivariate and the multivariate measures 
depend not only on the relative mixture of clusters but 
also on the total number of neurons. This has important 
implications for the analysis of multi-unit recordings that 
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typically access only a rather small number of neurons, 
and for which the relative amount of neurons that do or 
do not code for the presented stimulus at that moment 
in time is not known. 

We used a network of simulated Hindemarsh-Rose neu- 
rons to carry out a controlled comparison between the 
ISI- measures and six previously published methods. In 
this comparison the averaged bivariate measures perform 
better than the multivariate measures in distinguishing 
different sets of spike trains. Although the averaged bi- 
variate ISI-distance is not the overall best performer (this 
is the averaged bivariate Victor-distance), it is the mea- 
sure which combines sensitivity and time-resolved visu- 
alizability best. Among the multivariate measures, the 
multivariate ISI-diversity is indeed the best performer. 
This is most likely due to the insensitivity of the other 
multivariate measures, the Tiesinga and the Hunter et al. 
reliabilities, to the spike train of origin. Their first step 
is to pool the different spike trains into one cumulative 
spike train. Thus, for any single spike, the information 
about the spike train of origin is lost, i.e., the value of the 
reliability methods is invariant to shuffling spikes among 
the spike trains. Many different spike trains would yield 
the same pooled spike train (in principle even one spike 
train containing all of the spikes plus many empty spike 
trains) and thus much of the information distinguishing 
different sets of spike trains taken from different clusters 
is lost. 

This limitation does not affect the multivariate ISI- 
diversity, which is more sensitive to local structure within 
the single spike trains, since the interspike intervals are 
defined by neighboring spikes of the same spike train. 
However, the multivariate ISI-diversity still performs not 
as well as the averaged bivariate measures, because the 
calculation of the instantaneous mean and standard de- 
viation across spike trains implicitly assumes that a uni- 
modal distribution exists. Spike trains belonging to dif- 
ferent clusters result in a multimodal distribution which 
cannot be well-described by a common mean and a com- 
mon standard deviation. Thus, although the multivari- 
ate ISI-diversity is a step in the right direction, further 
improvements are needed for it and for all multivariate 
methods in order for them to compete with the averaged 
bivariate methods. 

One difference between the ISTmeasures and most of 
the other methods is that for the former, we neither have 
to divide the time series into windows of a certain fixed 
size, nor need to choose a parameter that sets the time 
scale of analysis. In our quantitative comparison this dif- 
ference was less apparent since we only showed results 
obtained for time scale parameters that have already 
been optimized. However, in applications to real data 
for which there is no validated knowledge about the data 
and its relevant time scales, it can be preferable to get a 
more objective estimate of neuronal variability by using 
a method that is time scale independent. Of course the 
computational cost and the time and effort that is needed 
to find the right parameter, if there indeed is one, might 



also play a role. However, it is also important to keep in 
mind that for studies dealing with specific questions of 
neuronal coding (e.g., what is the time scale over which 
different stimuli can be best distinguished?) the use of 
a parameter-based methods is preferable. But even in 
such cases, it might be useful to get a time scale inde- 
pendent estimate of the fundamental distinguishability of 
the spike trains in question. This holds particularly true 
for spike trains that contain different time scales such 
as regular spiking and bursting, since in these cases the 
optimum time scale can become ambiguous. 

In the third part of this study we reanalyzed single- unit 
monkey data to compare our methods with the Fano fac- 
tor, a standard method that employs a moving window 
technique. In this application we qualitatively demon- 
strated that our instantaneous measures allow us to an- 
alyze the data at the maximum possible resolution, but 
that this resolution can then easily be reduced by means 
of a moving average filter with the desired order. On the 
other hand, for measures like the Peri-Stimulus-Time- 
Histogram (PSTH) and the Fano factor, which rely on 
windowing, maximum resolution is always limited due to 
the demands on the spike counting statistics. 

A few of the caveats that one should take into account 
for the analysis of two spike trains also hold for the mul- 
tiple spike train case. First, it is obvious that no measure 
that results in a single number quantifying the synchrony 
between two or more spike trains can be adequate to 
deal with all kinds of potential coding schemes (e.g, time 
coding, rate coding and pattern coding; cf. also Vic- 
tor and Purpura [29|. Then, the ISTmeasures, like the 
other measures, are not sensitive to phase lags. Thus any 
such phase lags should be removed by suitably shifting 
the time series before computing the measures. Methods 
for the detection of lag synchronization in neuronal data 
[sot well as for their elimination [20] have been de- 
scribed. Finally, for the ISI-distances as for most of the 
other measures relative changes in value across different 
sets, stimuli, conditions, etc. are often more important 
than the absolute values. 

Note that the multiple spike train extensions proposed 
in this study aim at the quantification of synchrony 
within one set of spike trains. For the Victor-Purpura 
and the van Rossum distances, there also exists a second 
group of extensions designed to esti mate the synchr ony 
between two popul ations of neurons (jAronov 1 . 1200 3 and 
iHoughton and 2008, respectively). A correspond- 
ing extension of the ISI-distance will be presented in a 
forthcoming study [iGl] . 

In summary, due to their time-resolved defini- 
tion, time scale independence and computational 
speed, we expect the ISI-measures to be a powerful 
tool for all kinds of applications on multiple spike 
train synchrony and reliability. We close by point- 
ing out once more that the Matlab source codes 
for calculating and visualizing the instantaneous fir- 
ing rate and the ISI-measures can be downloaded at 
http://inls.ucsd. edu/^kreuz/Source-Code/Spike-Sync.html[ 
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APPENDIX A: ANALYTICAL RESULTS FOR 
POISSON PROCESSES 

In this Section we derive analytical results for the mean 
and the standard deviation for the averaged bivariate ISI- 
distance Dj and the ISI-diversity DJ^ in the case of mul- 
tiple Poisson processes. For this process the spike count 
in a window follows a Poisson distribution, whereas the 
ISI distribution is exponential. 

An exponential ISI-distribution as obtained for a Pois- 
son process with rate A yields a coefficient of variation 
Cv — 1 (both standard deviation and mean are equal to 
1/A). However, relevant for the instantaneous measures 
is the probability at each time instant to observe the in- 
terspike interval x. This probability is the product of the 
length and the frequency of the interval (with both of 
these quantities expressed in units of the mean 1/A), i.e., 



P{x) = X^xe 



(Al) 



X2e 



-\X2 
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Note that the result is independent of the rate. According 
to Eqs. [5] and [7] the same expectation value is obtained 
for the averaged bivariate ISI-distance Dj and the in- 
stantaneous average A{t). 

Correspondingly, for the ISI-diversity D™ we calculate 
the expectation value of the instantaneous coefficient of 
variation Cy(t) which we get from the mean and the 
standard deviation of the above distribution (Eq. 

The mean of this distribution is 



dxX^x^e 
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whereas its standard deviation is 
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This yields the expectation value for the coefficient of 
variation 



E{Cv{t)) = 



1 



(A5) 



which is also independent of the rate but in fact consid- 
erably smaller than 1. 

In Fig. [T7]we show the instantaneous analysis for 300 
Poisson spike trains each with the same rate A = 1 . The 
expected values < ISI >= 2 (Eq. EH with A = 1), 
A = 0.5 and Cv = l/\/2 are marked by dashed red lines. 
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For the bivariate ISI-distance Dj we calculate the expec- 
tation value of the instantaneous ISI-ratio I{t) (Eq. [3]) 
using this probability distribution for both xi and X2 [j]: 

E{I{t)) ^X^ / dxidx2{l - —)xie-^'''x2e-^''^' -fA2) 

Jo Jo ^2 

/ / dxidx2{l - —)xie-^''\ 
Jo Jo 
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FIG. 17: (color online) Multiple Poisson spike trains {N = 300) 
each with the same rate A = 1. Below are shown the instanta- 
neous mean ISI, the instantaneous average bivariate ISI-ratio A 
and the multivariate coefficient of variation Cv- The dashed red 
lines mark the expectation values for the instantaneous variability 
estimates calculated according to Eqs. IA2l and [A5l 



APPENDIX B: DATA 

All recordings and simulations were performed prior to 
and independently from the design of this study. 



1. In vitro recordings of a cortical cell 

In Section [ni] the ISI-measures are illustrated using 
in vitro whole-cell recordings taken from a cortical cell 
from the layer 2 of the medial entorhinal cortex of a 
young Long-Evans rat. In these experiments (conducted 
as approved by the University of California at San Diego 
Institutional Animal Care and Use Committee) cortical 
cells were selected from a 400 micron slice preparation 
by their superficial position, as well as particular char- 
acteristics of their electrophysiological responses to long 
current steps [cf. P^. Intracellular signals were ampli- 
fied, low pass filtered, and digitized at 10 kHz via soft- 
ware created in Lab View (National Instruments, Austin, 
TX, USA). Inputs were delivered as synaptic conduc- 
tances through a linux-based dynamic clamp Inputs 
were comprised of synaptic inputs added to an under- 
lying DC depolarization. The amplitude of the DC de- 
polarization was adjusted for each cell to elicit a spike 
rate of 5 — 10 Hz. Synaptic inputs were of the form 
Ixnn = GsynSiVm — Vsyn) where S followed the differ- 



S) - pS; a ^ 500/ms, 
was tailored for each 



ential expression dS/dt ~ a{\ 
f3 = 250/ms, Vsyn = mV. Ggyn 
cell to be peri-threshold and elicit a spike with probability 
close to 50%. Synaptic events were delivered for 10 sec- 
onds at irregular intervals that were taken from record- 
ings of spike times in response to steady DC depolariza- 
tion alone. Inputs were temporally scaled resulting in 
firing rates that were roughly matched to the elicited fir- 
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ing rate of the neuron, about one third slower and about 
twice as fast, with the latter covering only about half 
the length of the recording time. For each scaling five 
repetitions were analyzed. 

2. Hindemarsh-Rose simulations 

The spike trains of the controlled configuration used 
in Section [IVI were generated using time series extracted 
from a network of Hindemarsh-Rose neurons [9i]. This 
network was originally designed to analyze semantic 
memory representations using feature-based models. De- 
tails of the network architecture and the implementation 
of the feature coding can be found in Morelli et al. [l^ . 
The clusterii ig p roperties of these data were detailed in 
Kreuz et al. [3- 

In short, the state of the neuron n was determined 
by three dimensionless first-order differential equations 
describing the evolution of the membrane potential X„, 
the recovery variable y„, and a slow adaptation current 

Xn - Yr,-X^+3X^^-Z„+In+an-Pn (Bl) 

y„ = l~5X^-Y„ (B2) 
Z„ = 0.006[4(X„-1.6)-Z„], (B3) 

where 

F([/-l) 

a„ = ^ WnmArn (B4) 
m— 1 

and 

1 

k=l 

are the weighted inter-modular and intra-modular activ- 
ities, respectively, which were dependent on the synaptic 
connection weights Wnm and the neuronal activity vari- 
ables An. The external input /„ was chosen such that 
the Hindemarsh-Rose neurons were operating in a chaotic 
regime. 

The network consisted of 128 Hindemarsh-Rose neu- 
rons belonging to C/ = 16 different modules with F = 8 
neurons each. In a learning stage, input memory pat- 
terns were stored by updating the synaptic connection 
weights Wnm between different neurons using a Hebbian 
mechanism based on the activity variables An ■ A neuron 
was considered active whenever its membrane potential 
X exceeded the threshold value X = and its activity is 
coded by the variable An = 9{Xn — X), where 9 is the 
Heavyside function with 9{x) = 1 if x > and 9{x) = 
if a; < 0. For this study we restricted ourselves to 26 time 
series Xn extracted during the retrieval stage in which the 
learned connection weights were kept constant. Accord- 
ing to their coding properties regarding the retrieval of 



only two distinguished memory patterns, they belonged 
to two principal clusters: 13 of the neurons coded for 
pattern 1 only, and 13 coded for pattern 2 only. The 
respective time series were labelled by "1" and "2" fol- 
lowed by an index letter. The numerical integration was 
done using a fourth-order Runge-Kutta integration with 
a fixed step-size of 0.05 (arbitrary time units). The length 
of the time series was 32768 data points. The threshold 
for spike detection was chosen as the arithmetic average 
over the minimum and maximum value of the respective 
time series. 



3. Monkey data from a visual attention task 

In Section|V]our instantaneous measures are compared 
against standard methods of rate and variability estima- 
tion. This comparison is carried out on single-unit data 
from a neuron in area V4 of a rhesus macaque mon- 
key that were recorded during an attention-demanding 
multiple-object tracking task. These data were described 
in more detail in Mitchell et al. [l^. In short, at the be- 
ginning of the recording session, the receptive field of the 
neuron was mapped using subspace reverse correlation on 
60 Hz Gabor stimuli. The monkey began each trial by 
fixating a central point for 200 ms and then maintained 
fixation through the trial. Four identical Gabor stimuli 
appeared outside the neurons receptive field at equally 
eccentric positions separated by 90°. Two stimuli were 
briefly identified as targets before all stimuli moved along 
independent trajectories at approximately 10° /s for 950 
ms, placing them at a new set of equally spaced loca- 
tions, one of which was within the receptive field. The 
stimuli paused for 1000 ms before moving to another set 
of locations and stopping. In the end of the trial the 
fixation point disappeared, and the monkey made a sac- 
cade to each target. Trials were categorized as either 
attended or unattended depending on the type of stim- 
ulus that entered and paused in the receptive field. For 
each of the two attention conditions eight novel sets of 
stimulus trajectories were generated and each repeated 
six times resulting in a total number of 96 trials. The 
data can be downloaded from J. F. Mitchell's webpage: 
,http://www.snl.salk.edu/^jude/sfn2008/index.html, 



APPENDIX C: PREVIOUSLY PUBLISHED 
MEASURES OF SPIKE TRAIN DISTANCE 

In Section IIVI we use a controlled configuration to 
compare the performance of the average bivariate ISI- 
distance Dj and the multivariate ISI-diversity D™ to four 
averaged bivariate and two multivariate measures. The 
four averaged bivariate measures Dy, Z)^, Dg and Dq 
are based on the spike train distances introduced by Vic- 
tor and Purpura [28[ and van Rossum [23, the similar- 
ity measure proposed by Schreiber et al. [22l |. and event 
synchronization [21 ] , respectively. Multivariate measures 
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comprise the reliabilities by Hunter et al. [T2| and 
Dip by Tiesinga 

The four measures Dy, _D^, Dg and rely on a pa- 
rameter that sets the time scale of the analysis. Following 
Kreuz et al. these parameters were varied over several 
logarithmic decades, with four equidistant values within 
each decade, and then optimized on the Hindemarsh- 
Rose data, this time with respect to the present task 
(see Section ITV]). For the averaged bivariate methods the 
time scales yielding the best performance were identical 
to those found for the bivariate measures in Kreuz et al. 
[l^ . where the performance in reproducing the cluster- 
ing of the same Hindemarsh-Rose network was evaluated. 
This is not surprising since both of these tasks rely on 
the distinction of the same spike trains. 



1. Averaged bivariate measures 

For this class of measures, to which the averaged bi- 
variate ISI-distance belongs as well, synchrony is defined 
as the average value over all pairs of spike trains. The 
underlying bivariate measures are described below. To 
simplify notation the two spike trains are here labelled x 
and y. 



Since the post-synaptic currents triggered by single spikes 
are well fit by exponentials, the van Rossum distance 
estimates the difference in the effect of the two trains 
on the respective synapses. In this method, the time 
constant of the exponential as the parameter that sets 
the time scale. Following optimization it is set to = 
177.83 time units. 



c. Schreiber et al. similarity measure 



In this corr elation-based approach ( Schreiber et al.l . 
120031 : cf. also iHaas and Whitd . l2002h each spike train 
is convolved with a Gaussian filter of width as to form 
x' and y' before cross correlation and normalization: 



C'sias) 



\x'm 



(C2) 



To derive a normalized measure of spike train dissimilar- 
ity the measure is inverted according to Ds = 1 — Cs- 
The width of the convolving filter ag sets the time scale 
of interaction between the two spike trains. Optimization 
resulted in a value of as = 31.62 time units. 



d. Event synchronization 



a. Victor- Purpura spike train distance 

The spike train distance Dy introduced in Victor 
and Purpura ^28*1 defines the distance between two spike 
trains in terms of the minimum cost of transforming one 
spike train into the other by using just three basic opera- 
tions: spike insertion, spike deletion and spike movement. 
While the cost of insertion or deletion of a spike is set to 
one, the cost of moving a spike by some interval is the 
only parameter of the method, and sets the time scale 
of the analysis. For zero cost, the distance is equal to 
the difference in spike counts, for high costs, the distance 
approaches the number of non-coincident spikes, as it be- 
comes more favorable to delete all non-coincident spikes 
than to shift them. Thus, by increasing the cost, the dis- 
tance is transformed from a rate distance to a timing dis- 
tance. Optimization of the cost on the Hindemarsh-Rose 
data yielded an intermediate value of 0.01/time unit. 



b. Van Rossum spike train distance 

A second spike train distance was introduced in van 
Rossum In this method, each spike is convolved with 
an exponential function with time constant r^. From 
the convolved waveforms x(t) and y{t), the van Rossum 
distance D_r can be calculated as 



1 r°° 

Dr{tr) = — / [x{t) - y{t)fdt 
tr Jo 



(CI) 



E vent synchronization ( Ouian Quiroga et al.l[2003 cf. 
also iKreuz et~ ], l2nn7bD works as a coincidence detec- 
tor quantifying the level of synchrony from the num- 
ber of quasi-simultaneous appearances of spikes. In con- 
trast to the measures introduced above, this method is 
parameter- and scale-free, since the time lag up to 
which two spikes tf and are considered to be syn- 
chronous is adapted to the local spike rates: 



r,-., = mm 



{^?Vl-^^t^-d,<?+l-<^^^^^l}/2. (C3) 



The function c{x\y) counts the appearances of a spike in 
X shortly after a spike in y: 



My 

1=1 j=l 



(C4) 



where 



if 0<if- 
if = 

otherwise. 



< T- 
3 — y 



(C5) 



With c{y\x) defined accordingly, event synchronization 
can be written as 



Q 



c{y\x) + c{x\y) 



(C6) 



A suitably normalized distance is obtained by inversion 
Dq = (5—1, with Dq = if and only if all spikes coincide. 



14 



Multivariate measures 



pooled spike train 



The first step for the Hunter et al. [I4I and the Tiesinga 
[25| reliability measures is to pool the spike times of 
all neurons together into one set i^, tf , t£^p with 
denoting the total number of spikes. The inter- 
spike intervals of this pooled spike train are labelled as 



tf 



a. Hunter et al. rehabihty 

Like the van Rossum distance the multivariate relia- 
bility measure introduced by Hunter et al. [H] takes 
distances between spikes into account by convolving each 
spike of the pooled spike train with an exponential func- 
tion with time constant th- Then it evaluates the vari- 
ance of the resulting continuous function, which is 
not bounded and expected to be high for synchronous 
and low for independent spike trains. The normalization 
proposed in Hunter et al. |l2| relies on some general as- 
sumptions that are not fulfilled in our case, e.g., the same 
number of spikes M in all spike trains. Thus, we normal- 
ized the variance by the value V,^^^ obtained for the set of 
identical spike trains that is generated by reproducing N 
times the single spike train of the set that gives the max- 
imum variance value, i.e., V^^^ = max„=i^...^Ar{y"} 
where denotes the variance of the train from neuron 
n. We then take the complement to 1 in order to get a 
measure of dissimilarity: 



o-(Pisi) 
(Pisi)^ ' 



(C8) 



For an asynchronous ensemble each neuron fires uncorre- 
lated to the other. The spike time histogram is therefore 
flat. For a homogeneous Poisson spike train, the coeffi- 
cient of variation Cy of a single neurons interspike inter- 
val distribution will be asymptotically equal to one. The 
superposition of independent Poisson processes is also a 
Poisson process with a Cy equal to one as well. The vari- 
ance in the Cy value across different realizations of the 
ensemble depends only on the number of spikes in 
the pooled spike train, not on the number of neurons N in 
the ensemble. By normalizing the Cy to the analytically 
computed value for a perfectly synchronous ensemble [for 
details see[25i] and inverting one obtains the measure of 
dissimilarity 



N 



(C9) 



that is sensitive to the precision of the ensemble discharge 
as well as to the degree of coincidence. This method is 
free of parameters. 



D 



H 



VP 



(C7) 



The time constant th of the exponential sets the time 
scale. Best results were obtained for th = 56.23 time 
units. 



b. Tiesinga reliability 

The first step for the measure introduced in Tiesinga 
I25I is to calculate the coefficient of variation for the 
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